function simu_cn0113a(lens)
% simulate copy number variant
% function simu_cn0113(lens)
%
% input:
% lens is the length of sequence
%
% output:


fid = fopen('1.simu_cn.names.txt','w');
fid1 = fopen('1.simu_cn.txt','w');


times=[3, 6, 9, 20, 50, 100];
times_len=length(times);
max_num=10000;

states=[0:4];
states1=[0,1,3,4];
states1=states1';
len_states1=length(states1);

lens=[2489565];

states_seq=2*ones(lens, 1);
rc=ones(lens, 1);

num_part=4;
max_part=.9;
max_part_1=.4;
min_part=.1;
min_part_1=.5;

zero=0;

for i=1:times_len
	% rands=randi(lens, times(i), 1);
	% [sigma00, sigma11, sigma22, sigma33, sigma44, sigma55]=ndgrid()
	sigma_states=.2*ones(1,1,len_states1);
	part04=min_part:(max_part_1-min_part)/num_part:max_part_1;
	part13=min_part_1:(max_part-min_part_1)/num_part:max_part;
	[part00, part11, part33, part44]=ndgrid(part04, part13, part13, part04);
	part_04=[part00(:), part11(:), part33(:), part44(:)];
	for j=1:length(part_04)
		part=part_04(j,:);
		obj=gmdistribution(states1, sigma_states, part);
		y=random(obj, times(i));
		y=round(y);
		y=y(find(y==0 |y ==1 | y==3 |y==4));
		len_y=length(y);
		rands=randi(lens, len_y, 1);
		% states_seq(rands(1:len_y))=y;
		
		stat1=randperm(6);
		lenx=3;
		stat1=stat1(1:lenx);
		stat1=stat1';
		sigma_stat1=.2*ones(1,1,lenx);
		[partt1, partt2, partt3]=ndgrid(min_part:(max_part-min_part)/num_part:max_part);
		partt=[partt1(:), partt2(:), partt3(:)];
		
		for k=1:length(partt)
			parttt=partt(k,:);
			obj1=gmdistribution(stat1, sigma_stat1, parttt);
			z=random(obj1, len_y);
			z=round(10.^z);
			for ii=1:len_y
				abnorm_start=rands(ii);
				abnorm_end=rands(ii)+z(ii);
				if abnorm_end > lens
					states_seq(abnorm_start:lens)=y(ii);
				else
					states_seq(abnorm_start:abnorm_end)=y(ii);					
				end
			end
			
			x0=find(states_seq==0);
			x1=find(states_seq==1);
			x2=find(states_seq==2);
			x3=find(states_seq==3);
			x4=find(states_seq==4);
			
			stat01=[0,1]';
			sigma_stat2=ones(1,1,2);
			part22=[.95, .05; .85, .15; .75, .25];
			for jj=1:size(part22)
				part2=part22(jj,:);
				obj2=gmdistribution(stat01, sigma_stat2, part2);
				if length(x0)>0
					rcx=random(obj2, length(x0));
					rc(x0)=rcx;
				end
				
				stat012=[0,1,2]';
				sigma_stat012=ones(1,1,3);
				part012=[.025, .95, .025; .1, .8, .1; .15, .7, .15];
				for jj=1:size(part012)
					part2=part012(jj,:);
					obj2=gmdistribution(stat012, sigma_stat012, part2);
					if length(x1)>0
						rcx=random(obj2, length(x1));
						rc(x1)=rcx;
					end
					
					stat123=[1,2,3]';
					for jj=1:size(part012)
						part2=part012(jj,:);
						obj2=gmdistribution(stat123, sigma_stat012, part2);
						if length(x2)>0
							rcx=random(obj2, length(x2));
							rc(x2)=rcx;
						end
						
						stat234=[2,3,4]';
						for jj=1:size(part012)
							part2=part012(jj,:);
							obj2=gmdistribution(stat234, sigma_stat012, part2);
							if length(x3)>0
								rcx=random(obj2, length(x3));
								rc(x3)=rcx;
							end
							
							stat43=[4,3]';
							for jj=1:size(part22)
								part2=part22(jj,:);
								obj2=gmdistribution(stat43, sigma_stat2, part2);
								if length(x3)>0
									rcx=random(obj2, length(x4));
									rc(x4)=rcx;
									rc(find(rc<0))=0;
									
									zero=zero+1;
									
									if max_num < lens
										rc1=randsample(rc, max_num);
										len1=max_num;
									else
										rc1=rc;
										len1=lens;
									end

									rc1=sort(rc1);
									sum1=sum(rc1);

									for kk=1:len1
										rc_norm(kk)=sum(rc1(1:kk))/sum1;
									end
									gini=1-((2*sum(rc_norm)-rc_norm(1)-rc_norm(len1))/(len1-1));
									
									name=['1.', '.', num2str(zero), '.simu_cn.txt'];
									fprintf(fid, '%s %.3f\n', name, gini);
									dlmwrite(name,rc);
								end
							end
						end
					end
				end
			end
		end
	end
end



